# 	Filename: exitLD_source.R
# 	Description: Source script to produce results based on regressions where `leaderexitsoffice' is the dependent variable using the data unimputed set and listwise deletion. Called from exit_driver.R
# 	Author: Barry Hashimoto
# 	Date: December 2019
# 	For: Barry Hashimoto, "Autocratic Consent to International Law: the Case of the International Criminal Court's Jurisdiction," International Organization.

############################################################################################################
# Prep data.
 data <- NULL
 data <- all
 data$ruleoflaw <- data$vdem.ruleoflaw
 library(caret)
 legaldummies <- as.data.frame(predict(dummyVars(~ as.factor(legaltrad2), data = data), newdata = data))
 names(legaldummies) <- c("civil", "common", "mixed", "islam")
 data <- data.frame(data, legaldummies)
 data$ccode <- factor(data$ccode)
 
# Assign all countries with max(leaderexitsoffice) = 0 the same country code for efficient Firth logit estimation.
  NewCcode <- with(data, aggregate(leaderexitsoffice, by = list(ccode=ccode), FUN = max))
  NewCcode$x <- as.factor(ifelse(NewCcode$x==0, 999, NewCcode$ccode))
  names(NewCcode) <- c("ccode", "ccode2")
  data <- merge(data, NewCcode, by = "ccode")
 
# Identify the OECD DAC States so they can be dropped
  dacmemberlist <- c("Australia", "Austria", "Belgium", "Canada", "Denmark", "Finland", "France", "Germany", 
    "Greece", "Ireland", "Italy", "Japan", "South Korea", "Luxembourg", "Netherlands", "New Zealand", "Norway", 
    "Portugal", "Spain", "Sweden", "Switzerland", "United Kingdom", "United States of America")  	    
 data <- data[ifelse(data$country %in% dacmemberlist, 1, 0)==0,]
 data <- data[data$quarter>=154, ]

 writeLines("")
 writeLines("Print the total N of observations for both autocracies and democracies since 1998 Q3, before any listwise deletion.")
 print(dim(data)[1])
 
############################################################################################################
# Created matched data sets. RegimeOption is set in the file exit.driver.R to equal 0 for autocracies, 1 for democracies.
  set <- subset(data, subset = democracy == RegimeOption)

 writeLines("")
 writeLines("Print the total N of observations for only the selected regime type since 1998 Q3, before any listwise deletion.")
 print(dim(set)[1])
   
# Country-specific time-trends estimated from cubic polynomial
  set$ccode <- factor(set$ccode)
  set$exit.trend.single <- predict(lm(leaderexitsoffice ~ poly(quartersinoffice,3), data = set))
  set$exit.trend.country <- 10*predict(lm(leaderexitsoffice ~ ccode*poly(quartersinoffice,3), data = set))

  data.list <- with(set, data.frame(leaderexitsoffice, logdeaths, ptssum.mean, RomeStatuteRatification, 
    logoilrent, allunemp, loggrowth, x.ucdp.incidence, capital.scaled, gdp.scaled, ruleoflaw, democracy, 
    common, mixed, islam, pre98conflict, medianIMR, ccode, year, exit.trend.single, exit.trend.country, lcode))
  
# Complete cases.
  data.list <- data.list[complete.cases(data.list), ]
  length(unique(data.list$ccode)); length(unique(data.list$lcode))

# Prep for matching. 
  breaks <- 1/3
  oil.cuts <- quantile(data.list$logoilrent, seq(0,1, breaks))
  allunemp.cuts <- quantile(data.list $allunemp, seq(0,1, breaks))
  growth.cuts <- quantile(data.list $loggrowth, seq(0,1, breaks))
  gdp.cuts <- quantile(data.list $gdp.scaled, seq(0,1, breaks))
  capital.cuts <- quantile(data.list $capital.scaled, seq(0,1, breaks))
  rol.cuts <- quantile(data.list $ruleoflaw, seq(0,1, breaks))
  pts.cuts<- quantile(data.list $ptssum.mean, seq(0,1, breaks))
  imr.cuts <- quantile(data.list $medianIMR, seq(0,1, breaks))

 cut.list <- list(allunemp = allunemp.cuts, loggrowth = growth.cuts, 
    logoilrent = oil.cuts, gdp.scaled = gdp.cuts, capital.scaled = capital.cuts, 
    ruleoflaw = rol.cuts, ptssum.mean = pts.cuts, medianIMR = imr.cuts)

# A function to stop excess verbosity when running matchit().
 quiet <- function(x) { 
   sink(tempfile()) 
   on.exit(sink()) 
   invisible(force(x)) 
  }
  
# Match data by CEM. Note: cowwar dropped due to no variance for either dems or auts.
 matched.data <- quiet(matchit(RomeStatuteRatification ~ x.ucdp.incidence + logoilrent 
     + allunemp + loggrowth + gdp.scaled + ruleoflaw + capital.scaled + common + mixed 
     + islam + pre98conflict + medianIMR + ptssum.mean, data = data.list, method = "cem", 
     cutpoints = cut.list))
     
############################################################################################################
# Full data, logit, shared baseline, slimmed model without controls
 shared.nocontrol <- glm(formula = leaderexitsoffice ~ RomeStatuteRatification + exit.trend.single
   + as.factor(ccode2) + as.factor(year), family = binomial(link = "logit"), data = set)

# Full data, logit, shared baseline with controls
 shared.control <- glm(formula = leaderexitsoffice ~ RomeStatuteRatification 
    + x.ucdp.incidence + logoilrent + loggrowth + allunemp + gdp.scaled + ruleoflaw 
    + capital.scaled + exit.trend.single + as.factor(ccode2) + as.factor(year), 
    family = binomial(link = "logit"), data = set)
	
# Matched data, logit, shared baseline
 shared.matched <- glm(formula = leaderexitsoffice ~ RomeStatuteRatification 
   + x.ucdp.incidence + logoilrent + loggrowth + allunemp + gdp.scaled + ruleoflaw 
   + capital.scaled + common + mixed + islam + pre98conflict + medianIMR + ptssum.mean 
   + exit.trend.single, family = binomial(link = "logit"), data = match.data(matched.data))

# Full data, logit, unit-specific baseline, slimmed model without controls
  country.nocontrol <- glm(formula = leaderexitsoffice ~ RomeStatuteRatification 
  + exit.trend.country + as.factor(ccode2) + as.factor(year), family = binomial(link = "logit"), 
  data = set)

# Full data, logit, unit-specific baseline with controls 
  country.control <- glm(formula = leaderexitsoffice ~ RomeStatuteRatification + x.ucdp.incidence 
     + logoilrent + loggrowth + allunemp + gdp.scaled + ruleoflaw + capital.scaled + exit.trend.country 
     + as.factor(ccode2) + as.factor(year), family = binomial(link = "logit"), data = set)

# Matched data, logit, unit-specific baseline
  country.matched <- glm(formula = leaderexitsoffice ~ RomeStatuteRatification + x.ucdp.incidence 
    + logoilrent + loggrowth + allunemp + gdp.scaled + ruleoflaw + capital.scaled + common + mixed 
    + islam + pre98conflict + medianIMR + ptssum.mean + exit.trend.country, 
    family = binomial(link = "logit"), data = match.data(matched.data))

# Functions for calculating estimates which works with glm(...logit).
 simit <- function(x){
  require(mvtnorm)
  betas <- summary(x)$coefficients[,1]
  variances <- as.matrix(vcov(x))
  sim.mat <- rmvnorm(100000, betas, variances)
  sim.vec <- sim.mat[,2]
  estimates <- c(mean(sim.vec), sd(sim.vec), dnorm(mean(sim.vec)/sd(sim.vec)))
  invisible(estimates)
  }
 
 output <- function(x){
	mat <- rbind(simit(x), 
	             c(AIC(logLik(x)), NA, NA),
	             c(length(x$residuals), NA, NA)
	             )
	#mat[c(1:2),] <- signif(mat[c(1:2),], 4)               
	rownames(mat) <- c("Coef. of ICC Jurisdiction", "Model AIC", "N")
	return(mat)
	}

# SATT estimate.

 satt.exit <- zelig(formula = leaderexitsoffice ~ RomeStatuteRatification + x.ucdp.incidence + logoilrent + loggrowth + allunemp + gdp.scaled + ruleoflaw + capital.scaled + common + mixed + islam + pre98conflict + medianIMR + ptssum.mean + exit.trend.country, model = "logit", data = match.data(matched.data), cite = F) %>% ATT(treatment = "RomeStatuteRatification", treated = 1, num = 20000) %>% get_qi(qi="ATT", xvalue = "TE")

 writeLines("")
 writeLines("Print the 95% interval plus 25th and 7th percentiles (top) and median absolute deviation (bottom) of the SATT of RomeStatuteRatification on leader exit calculated from models fit to the matched data w/listwise deletion for the selected regime type.")
 print(signif(quantile(satt.exit, c(0.025, 0.25, 0.5, 0.75, 0.975)), 2))
 print(signif(mad(satt.exit), 2))

   
######################################################################################################
# Table of results for `leaderexitsoffice'. Unimputed data.

 col1 <- rbind(output(shared.nocontrol), output(shared.control), output(shared.matched))
 col2 <- rbind(output(country.nocontrol), output(country.control), output(country.matched))
 exitable <- cbind(col1, col2)
 
 colnames(exitable) <- rep(c("Est.", "SE", "p"), 2)

 writeLines("")
 writeLines("Print the table of results for models of leader exit using the unimputed data for the selected regime type, with listwise deletion. First column has models with shared trends, while second column has models with country-specific trends. First row has results for unmatched data with no control variables. Second row has results for unmatched data with control variables and the same fixed effects. Third row has results for models using the matched data sets with control variables.")
 print(round(exitable, 3))
 dim(data.list)
 
 
 
 
 
 
 
